#-*- coding: utf8 -*-
from tutor.scripts import *
from tutor.plugins.maple import *

# variáveis úteis
half = M('1/2')
t = M.t
_i, _j, _k = M('_i, _j, _k')
rr = M('`..`')
coeff = M.coeff
def sumt(a, b): 
    return tuple( ai + bi for (ai,bi) in zip(a, b) )
def mul(c, v): 
    return tuple( c * vi for vi in v )
u, v = M.u, M.v

# parâmetros iniciais - vetores
a, b = M(oneof(1, 1, 2, 3)), M(oneof(1, 1, 2, 3))
Pk1 = (a, b, 0)
a, b = M(oneof(1, 1, 2, 3)), M(oneof(1, 1, 2, 3))
Pk2 = (a, 0, b)
a, b = M(oneof(1, 1, 2, 3)), M(oneof(1, 1, 2, 3))
Pk3 = (0, a, b)
a = M(oneof(1, 1, 2, 3))
Pk4 = (0, 0, a)

# vetor de campo
dir = oneof(_i, _j, _k)
var = oneof(x, y, z)
F = fieldF = var * dir
Fx, Fy, Fz = coeff(F, _i), coeff(F, _j), coeff(F, _k) 
print('F: %s' % F)

# vértices 
Pv1 = Pk1
Pv2 = tuple( ki+vi for (ki, vi) in zip(Pk1, Pk2) )
Pv3 = tuple( ki+vi for (ki, vi) in zip(Pk1, Pk3) )
Pv4 = tuple( ki+vi+ri+si for (ki,vi,ri,si) in zip(Pk1, Pk2, Pk3, Pk4) )
print('k1: %s\nk2: %s\nk3: %s\nk4: %s' % (Pk1, Pk2, Pk3, Pk4))
print('v1: %s\nv2: %s\nv3: %s\nV4: %s' % (Pv1, Pv2, Pv3, Pv4))

# cálculo do fluxo
X, Y, Z = sumt(sumt(sumt(Pk1,
                         mul(1-u, Pk2)),
                         mul(1-v, Pk3)),
                         mul((1-u)*(1-v), Pk4))

# derivadas da parametrizacao
Ru = M.Vector([diff(X, u), diff(Y, u), diff(Z, u)])
Rv = M.Vector([diff(X, v), diff(Y, v), diff(Z, v)])

# elemento de superficie
dS = M('LinearAlgebra:-CrossProduct')(Ru, Rv)
if subs(u==0.5, v==0.5, dS[3]) < 0:
    dS = -dS
print('dS: %s' % dS)

# integrando
JS = subs(x==X, y==Y, z==Z, dS[1] * Fx + dS[2] * Fy + dS[3] * Fz)

# integral de superfície
aux = M.int(JS, u==r(0,1))
Phi = M.int(aux, v==r(0,1))
print('fluxo:', Phi)

# erro: associar r2, r3, r4 a i, j, k
Xe, Ye, Ze = 1-u, 1-v, (1-u) * (1-v)
Rue = M.Vector([diff(Xe, u), diff(Ye, u), diff(Ze, u)])
Rve = M.Vector([diff(Xe, v), diff(Ye, v), diff(Ze, v)])
dSe = M('LinearAlgebra:-CrossProduct')(Rue, Rve)
if subs(u==0.5, v==0.5, dSe[3]) < 0:
    dSe = -dSe
JSe = subs(x==Xe, y==Ye, z==Ze, dSe[1] * Fx + dSe[2] * Fy + dSe[3] * Fz)
aux = M.int(JSe, u==r(0,1))
PhiE = M.int(aux, v==r(0,1))
print('fluxo errado:', PhiE)

# erro: confundir produto vetorial com escalar
dSe = M.Vector([Ru[1]*Rv[1], Ru[2]*Rv[2], Ru[3]*Rv[3]])
if subs(u==0.5, v==0.5, dSe[3]) < 0:
    dSe = -dSe
JSe = subs(x==X, y==Y, z==Z, dSe[1] * Fx + dSe[2] * Fy + dSe[3] * Fz)
aux = M.int(JSe, u==r(0,1))
PhiEsc = M.int(aux, v==r(0,1))
print('fluxo errado:', PhiEsc)

# erro: mistura os dois erros anteriores
dSe = M.Vector([Rue[1]*Rve[1], Rue[2]*Rve[2], Rue[3]*Rve[3]])
if subs(u==0.5, v==0.5, dSe[3]) < 0:
    dSe = -dSe
JSe = subs(x==Xe, y==Ye, z==Ze, dSe[1] * Fx + dSe[2] * Fy + dSe[3] * Fz)
aux = M.int(JSe, u==r(0,1))
PhiE3 = M.int(aux, v==r(0,1))
print('fluxo errado:', PhiE3)

PhiD = Phi + PhiE
print('fluxo errado:', PhiD)
